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Abstract 

We  consider  the  linear  stability  of  a modulated  Taylor-Couette  system  when 
the  inner  cylindrical  boundary  consists  of  a crystalline  solid-liquid  interface.  Both 
experimentally  and  in  numerical  calculations  it  is  found  that  the  two-phase  system 
is  significantly  less  stable  than  the  analogous  rigid-walled  system  for  materials 
with  moderately  large  Prandtl  numbers.  A numerical  treatment  based  on  Floquet 
theory  is  described,  which  gives  results  that  are  in  good  agreement  with  preliminary 
experimental  findings.  In  addition,  this  instability  is  further  examined  by  carrying 
out  a formal  asymptotic  expansion  of  the  solution  in  the  limit  of  large  Prandtl 
number.  In  this  limit  the  Floquet  analysis  is  considerably  simplified,  and  the  linear 
stability  of  the  modulated  system  can  be  determined  to  leading  order  through  a 
conventional  stability  analysis,  without  recourse  to  Floquet  theory.  The  resulting 
simplified  problem  is  then  studied  for  both  the  narrow  gap  geometry  and  for  the 
case  of  a finite  gap.  It  is  surprising  that  the  determination  of  the  linear  stability  of 
the  two-phase  system  is  considerably  simpler  than  that  of  the  rigid-walled  system, 
despite  the  complications  introduced  by  the  presence  of  the  crystal-melt  interface. 
PACS:  47.30. +s,47.20.-k,  81.10.Fq,81.30.Fb 
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1.  Introduction 

An  important  consideration  in  the  growth  of  multicomponent  single  crystals  from  the 
melt  is  the  spatial  distribution  of  solute.  In  many  instances,  inhomogeneities  in  the 
solute  distribution  result  in  inferior  electrical  and  mechanical  properties.  It  is  therefore 
desirable  to  understand  the  physical  mechanisms  that  lead  to  solute  segregation  in  melt- 
grown  crystals. 

Since  solid  state  diffusion  normally  occurs  at  rates  several  orders  of  magnitude  lower 
than  in  the  melt  phase,  the  distribution  of  solute  in  the  crystal  is  primarily  determined 
by  the  solute  profile  in  the  liquid  immediately  ahead  of  the  advancing  crystal-melt  inter- 
face; as  the  liquid  transforms  into  solid,  little  rearrangement  of  the  solute  distribution 
takes  place.  Thus  the  prediction  of  the  solute  distribution  in  the  crystal  requires  the 
understanding  of  solute  transport  in  the  melt  adjacent  to  the  crystal-melt  interface  [1,  2]; 
in  particular,  mechanisms  that  produce  lateral  segregation  of  solute  along  the  interface 
play  a major  role. 

During  directional  solidification  (see,  e.g.,  [3])  under  ideal  processing  conditions,  a pla- 
nar crystal-melt  interface  advances  at  constant  velocity  into  a quiescent  melt,  in  which 
case  the  crystal  that  is  produced  is  spatially  homogeneous.  Loss  of  ideality  can  arise  in 
many  ways;  two  of  the  more  extensively  studied  mechanisms  are  morphological  instabil- 
ity of  the  crystal- melt  interface  [4,  5],  which  causes  a planar  crystal- melt  to  develop  a 
corrugated  lateral  structure  with  accompanying  solute  redistribution,  and  hydrodynamic 
instabilities  in  the  melt  [6,  7,  8],  in  which  solute  inhomogeneities  are  associated  with  the 
features  of  the  secondary  flow  fields. 

Although  these  two  mechanisms  can  be  studied  separately,  it  is  important  to  un- 
derstand how  the  two  are  related.  Hydrodynamic  flows  occur  frequently  during  crystal 
growth  from  the  melt,  either  naturally  or  by  design.  It  is  desirable  to  be  able  to  predict 
whether  a particular  flow  will  generate  a strong  coupling  with  the  crystal-melt  interface, 
possibly  enhancing  morphological  instability.  It  is  similarly  important  to  know  whether 
the  presence  of  the  crystal-melt  interface  will  cause  significant  changes  in  the  hydrody- 
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namic  instabilities  that  occur  in  classical,  fixed-boundary  systems.  If  the  flow/interface 
interaction  is  particularly  strong,  of  course,  then  the  classification  of  the  underlying  insta- 
bility as  being  either  morphological  or  hydrodynamic  in  nature  may  be  overly  simplistic. 

A number  of  authors  have  considered  the  morphological  stability  of  the  crystal-melt 
interface  in  the  presence  of  various  types  of  flow  fields.  Examples  include  plane  Couette 
flow  [9,  10],  thermosolutal  convection  [11,  12],  plane  stagnation  flow  [13,  14],  rotating  disk 
flow  [15],  and  the  asymptotic  suction  profile  [16].  The  effect  of  a crystal- melt  interface 
on  the  hydrodynamic  stability  of  the  melt  has  also  been  considered  for  a variety  of 
flows,  including  Rayleigh-Benard  convection  [17],  thermosolutal  convection  [7,  18,  19], 
plane  Poiseuille  flow  [20],  the  asymptotic  suction  profile  [16],  thermally-driven  flow  in  an 
annulus  [21,  20],  and  steady  Taylor- Couette  flow  [22,  23]. 

Particularly  strong  flow/interface  interactions  were  observed  in  two  of  these  cases. 
For  flow  in  a vertical  annulus  with  radial  heating,  with  the  outer  annular  boundary 
consisting  of  a cylindrical  crystal- melt  interface,  Fang  et  al.  [21,  20]  found  that  the 
two-phase  system  is  destabilized  by  an  order  of  magnitude  compared  to  the  rigid- walled 
system.  The  material  studied  was  succinonitrile  (SCN),  which  has  a moderately  large 
Prandtl  number  of  about  23.  In  the  rigid-walled  system,  the  most  dangerous  mode 
under  these  conditions  is  an  axisymmetric  buoyant  instability.  The  most  dangerous 
mode  for  the  two-phase  system  is  given  by  a non- axisymmetric,  helical  instability  which 
ensues  at  much  lower  radial  temperature  differences;  moreover,  this  mode  appears  as  a 
destabilization  of  a shear  mode  which,  in  the  rigid-walled  system,  is  considerably  more 
stable  than  the  buoyant  mode.  Numerical  parameter  studies  of  the  linear  stability  of  the 
two-phase  system  show  that  the  helical  instability  is  dominant  for  crystals  with  melts 
that  have  large  Prandtl  numbers.  In  steady  Taylor- Couette  flow  [24],  another  strong 
coupling  occurs  if  one  of  the  cylindrical  boundaries  again  consists  of  the  crystalline  phase 
of  the  melt  contained  in  the  annular  region  [22,  23];  this  effect  is  also  found  to  provide 
significant  destabilization  of  the  flow  for  melts  with  high  Prandtl  numbers. 

In  this  paper  we  consider  an  unsteady  version  of  the  Taylor- Couette  experiment,  in 
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which  strong  destabilization  of  the  hydrodynamic  instability  is  observed  for  the  two- 
phase  system.  We  consider  the  linear  stability  of  a Taylor- Couette  geometry  in  which 
the  entire  system  undergoes  a torsional  oscillation  about  the  cylindrical  axis.  Preliminary 
experiments  for  this  problem  have  been  performed  using  SCN,  and  show  a marked  desta- 
bilization of  the  system  over  a range  of  frequencies  [25].  Numerical  studies  have  also  been 
performed  using  Floquet  theory  to  describe  the  linear  stability  of  the  system;  the  numer- 
ical results  show  similar  behavior  with  qualitative  agreement  between  the  theory  and  the 
preliminary  experiments.  Since  the  strong  flow/interface  coupling  is  found  to  be  most 
significant  for  moderately  large  Prandtl  numbers,  we  examine  here  the  asymptotic  limit 
Pr  1 . Somewhat  surprisingly,  we  find  that  in  this  limit  we  are  able  to  obtain  a quite 
tractable  description  of  the  linear  stability  of  the  two-phases  system  without  recourse 
to  Floquet  theory.  The  strong  destabilization  provided  by  the  crystal-melt  interface  for 
Pr  1 actually  results  in  a much  simpler  analysis  for  the  two-phase  system  than  that 
obtained  for  the  single-phase  system  with  rigid  walls  [26,  27]. 

In  the  following  section  we  describe  briefly  preliminary  experiments  which  have  been 
performed  for  this  problem.  This  is  followed  by  a description  of  a numerical  treatment 
of  the  linear  stability  problem  based  on  a pseudospectral  spatial  discretization  of  the 
linearized  governing  equations  with  a Floquet  analysis  in  the  time  variable.  The  approach 
is  described  both  in  the  simpler  narrow  gap  limit  of  the  governing  equations  [28],  and  for 
the  more  general  case  of  a finite  gap.  The  asymptotic  analysis  for  large  Prandtl  number 
is  described  for  both  the  narrow  and  finite  gap,  and  the  numerical  and  asymptotic  results 
are  compared  with  the  experimental  results.  Some  remarks  concerning  the  relevance  of 
the  findings  to  the  more  general  problem  of  predicting  the  strength  of  flow/interface 
interaction  are  given  in  the  conclusion. 

2.  Modulated  Taylor- Couette  Flow 

The  experiment  is  performed  using  succinonitrile;  the  relevant  material  properties  are 
given  in  Ref.  [21].  The  apparatus  consists  of  two  coaxial  cylinders  which  are  sealed 
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together  at  both  ends  to  provide  a single,  rigid  unit;  the  succinonitrile  occupies  the 
space  between  the  cylinders.  A radial  temperature  difference  is  maintained  across  the 
sample  by  circulating  a cooling  fluid  with  temperature  T0  inside  the  inner  cylinder,  while 
the  whole  apparatus  is  situated  in  a bath  held  at  a uniform  temperature  T2  > T0.  If 
these  temperatures  are  chosen  on  either  side  of  the  melting  point  Tm  of  SCN,  so  that 
To  < TM  < T2,  then  a cylindrical  crystal-melt  interface  forms  at  an  intermediate  radius 
R1.  A schematic  diagram  with  a corrugated  crystal- melt  interface  is  shown  in  Figure  1. 
To  produce  a modulated  Taylor-Couette  flow,  the  entire  apparatus  is  oscillated  about  the 
cylindrical  axis  of  the  system.  Accurate  control  of  the  oscillation  was  devised  by  using  a 
mechanical  linkage  to  the  servo-motor  of  a chart  recorder  to  drive  the  system,  with  the 
input  to  the  chart  recorder  provided  by  a signal  generator  to  allow  a wide  variation  in 
the  waveform  and  frequency  of  the  modulation. 

The  outer  diameter  of  the  smaller  cylinder  is  1.2  cm,  and  the  inner  diameter  of  the 
larger  cylinder  is  3.3  cm.  The  length  of  the  straight  portion  of  the  apparatus  is  24.1  cm, 
with  a further  2 cm  of  tapering  on  each  end  where  the  walls  are  joined  together.  The 
choice  to  perform  the  rigid-body  rotation  of  the  coaxial  cylinders  without  differential 
rotation  was  made  to  avoid  complications  with  seals  that  allow  slip  at  the  endwalls. 
The  design  of  the  drive  allows  variation  in  both  the  frequency  and  amplitude  of  the 
torsional  oscillation  of  the  system.  In  the  preliminary  experiments  performed  thus  far,  the 
stability  boundaries  were  determined  by  fixing  the  frequency  and  varying  the  amplitude 
of  the  oscillation,  using  a pure  sinusoidal  waveform  from  the  signal  generator.  Oscillation 
frequencies  up  to  a few  Hertz  can  be  produced  reliably  with  this  design. 

In  a typical  experimental  run,  a crystal-melt  interface  is  established  within  the  glass- 
walled  annulus  by  maintaining  a temperature  difference  of  a few  degrees  Kelvin  across 
the  gap,  with  the  mean  temperature  near  the  melting  point.  For  small  amplitudes  of 
modulation,  the  interface  is  cylindrical  in  shape  except  near  the  tapered  endwalls  of 
the  container.  At  larger  amplitudes,  the  interface  is  observed  to  develop  well-defined 
oscillations  along  the  axial  direction.  The  resulting  deformation  of  the  interface  appears 
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axisymmetric,  and  the  amplitude  of  the  interface  deformation  appears  steady  in  time.  By 
performing  a series  of  runs,  stability  limits  for  the  modulation  amplitude  associated  with 
the  breakdown  of  the  cylindrical  crystal-melt  interface  could  be  determined  as  a function 
of  the  oscillation  frequency.  As  described  below,  the  two-phase  system  is  found  to  become 
unstable  at  amplitudes  that  are  an  order  of  magnitude  lower  than  those  required  for  the 
instability  of  the  analogous  single-phase,  rigid-walled  coaxial  system. 

Rudimentary  flow  visualization  is  possible  by  using  a video  camera  to  provide  images 
of  the  system  from  a side  view.  By  examining  the  motion  of  a particle  in  the  melt,  it  is 
possible  to  determine  the  nature  of  the  flow  field  associated  with  the  instability  of  the 
crystal-melt  interface.  At  a given  instant  of  time,  the  flow  consists  of  Taylor  vortices 
stacked  al&ig  the  apparatus  in  the  axial  direction,  with  the  wavelength  of  the  interface 
deformation  given  by  the  extent  of  a pair  of  the  counter-rotating  cells.  In  Figure  2 
we  show  a multiple  exposure  photograph  that  shows  the  motion  of  an  irregularly-shaped 
particle  in  a single  cell.  The  multiple  exposure  was  taken  in  two  steps.  First,  a continuous 
recording  of  the  video  tape  was  edited  to  give  a sequence  of  individual  frames  that  are 
each  separated  in  time  by  one  period  of  the  oscillation  frequency.  In  the  resulting  film 
the  particle  is  observed  to  circulate  along  a closed  path  in  the  same  azimuthal  plane.  The 
multiple  exposure  photo  was  then  assembled  from  several  of  these  frames  to  illustrate  the 
path  of  the  particle;  the  positions  of  the  particle  in  Figure  2,  however,  are  not  sequential, 
and  were  chosen  to  give  a rough  indication  of  the  particle  trajectory. 

The  particle  moves  radially  inward  from  the  hot  outer  cylinder  towards  the  minimum 
in  the  interface  deformation  (largest  liquid  gap),  where  hot  fluid  flows  towards  the  crystal- 
melt  interface,  and  returns  toward  the  outer  cylinder  near  the  maximum  in  the  interface 
deformation.  The  deformation  is  consistent  with  the  average  temperature  of  the  fluid, 
since  the  heating  provided  by  the  outer  cylinder  relative  to  the  relatively-cooler  crystal- 
melt  interface  implies  that  the  interface  melts  back  under  the  impingement  of  warmer 
fluid,  and  bulges  out  under  the  influence  of  the  cooler  fluid.  The  wavelength  of  the  Taylor 
vortices  is  significantly  longer  than  that  observed  for  the  single-phase  system;  quantitative 
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comparisons  will  be  provided  below,  together  with  the  discussion  of  theoretical  results. 

3.  Theory 

The  determination  of  linear  stability  for  rigid-walled  systems  undergoing  modulated  os- 
cillations has  been  performed  by  a number  of  authors,  including  Riley  and  Lawrence  [29], 
Carmi  and  Tustaniwskyj  [30],  Barenghi  and  Jones  [31],  Kuhlmann  et  al.  [32],  Wu  and 
Swift  [33],  and  Murray  et  al.  [34].  Here  we  consider  the  pure  torsional  oscillation  of  a 
Taylor- Couette  apparatus  for  the  case  of  a two-phase  system  with  a crystal- melt  interface 
providing  one  of  the  bounding  surfaces  of  the  flow.  The  analysis  assumes  an  axisymmet- 
ric  geometry,  which  is  consistent  with  the  experimental  observations  for  the  two-phase 
system.  Both  narrow  and  wide  gap  formulations  are  considered.  For  simplicity,  gravity- 
induced  buoyancy  effects  and  similar  effects  related  to  the  radial  density  stratification 
are  ignored;  these  effects  were  considered  for  the  unmodulated  case  in  Ref.  [35],  and  were 
found  to  provide  only  minor  corrections  to  the  stability  results. 

3.1.  Wide  Gap 

The  rigid-body  torsional  oscillation  for  sufficiently  large  angular  speeds  will  produce  a 
centrifugal  instability  in  the  base  flow  [26,  27].  As  depicted  in  Figure  1,  we  employ  a 
cylindrical  coordinate  system  (r1 ,(/)',  z'),  and  consider  the  linear  stability  of  an  idealized 
base  state  that  extends  to  infinity  in  the  axial  direction  (primed  quantities  will  be  dimen- 
sional). The  rigid  inner  wall  is  at  r'  = Ro,  and  the  outer  rigid  wall  is  at  r'  = R2  . The 
axisymmetric  crystal- melt  interface  is  given  by  r ' — R\{z' , t')\  in  the  base  state  the  inter- 
face is  cylindrical,  with  a constant  radius  Ri(z',t')  = R^.  The  gap  Ri(z',t')  < r'  < R2 
contains  the  incompressible  melt  or  liquid  phase,  and  the  crystal  or  solid  phase  occupies 
the  remainder  of  the  domain.  The  entire  apparatus  is  oscillated  torsionally,  so  that  at 
the  undisturbed  crystal- melt  interface  the  azimuthal  velocity  is  given  by  v'  = 
and  at  the  outer  wall,  v'  = R2tti(t)]  here  f)i(t)  = Q cos(u/i')  and  is  the  amplitude 
of  the  oscillations.  The  radial  and  axial  velocities  at  the  melt  boundaries  are  given  by 
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u'  = w'  — 0,  respectively.  In  the  gap,  the  melt  velocity  is  divergence-free  and  satisfies 
the  Navier-Stokes  equations  in  cylindrical  coordinates.  We  assume  that  the  flow  is  ax- 
isymmetric,  as  observed  in  the  experiment.  We  scale  the  radial  and  axial  velocities  with 
v/L  where  L = R2  — R\  is  the  gap  width,  and  v is  the  kinematic  viscosity,  and  scale  the 
azimuthal  velocity  with  RiQ.  We  scale  the  spatial  coordinates  as  r'  = Lr , z'  — Lz,  scale 
the  time  as  w't'  = r,  and  scale  the  pressure  as  pv2 / L2 . The  deviations  of  the  temperature 
fields  in  the  melt  and  solid  from  the  equilibrium  value  established  by  the  crystal-melt 
interface  are  scaled  with  the  temperature  difference  AT  across  the  gap,  for  example, 


Tl 


T'l-Tm(1-TIRi) 

AT 


(i) 


where 


at  = t2  - tm{i  - r/Ri), 


(2) 


and  T is  a capillary  constant  related  to  the  depression  of  the  melting  point  by  interface 
curvature  (the  Gibbs-Thomson  effect).  The  crystal-melt  interface  (at  R in  the  base 
state)  is  located  at 

r'  = Rx  + Lh(z,t ) (3) 


when  disturbed.  In  dimensionless  form,  the  domain  is  r0  < r < r2,  where  r0  = 7/s/(l  —77) 
and  r2  = 1/(1  — 77),  with  77  = R1/R2  and  775  = R0/R2.  The  melt  occupies  the  region 
r\  -f  h(z,t)  < r < r2,  where  r\  = 77/(1  — 77).  In  the  linearized  theory,  the  disturbance 


of  the  crystal-melt  interface  location  h will  be  taken  close  to  zero,  and  the  interfacial 
conditions  will  be  referred  to  in  the  usual  way.  Neglecting  all  derivatives  in  the 
azimuthal  direction,  we  obtain  the  governing  equations  on  T\  < r < r2, 


du  du  du 

— 1-  u— — b vj— , 

dr  dr  dz  2(1—77) 


du  dw 
dr  dz 

’ Ta— 
r 


= 0, 


dv  dv  dv  uv 

CU—  + u—  + w—  H 

dr  dr  dz  r 

dw  dw  dw 
u—  + u—  + w— 
dr  dr  dz 


dp  ^ du  + l du  u ^ d2u 


dr 
dv 
dr2 


dr2  r dr 

1 dv  v d2v 
+ + 


dz 2’ 


dz 2’ 


r dr 

dp  dw  1 dw  d2w 
dz  dr2  ^ r dr  ^ dz2  ’ 


(4) 

(5) 

(6) 
(7) 
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where  the  Taylor  number  is  given  by  Ta  = 2 Ri$l2  L3  / v2 , and  uj  — uj'L2 jv\  here  v in  the 
kinematic  viscosity  of  the  melt.  In  addition  to  the  equations  of  fluid  motion,  there  will 
be  energy  equations  in  the  melt  and  solid.  In  the  melt, 

1 F)Tr  WTt  \ 

(8) 


otl  dTL  arL 

uj  — b U— b W- 


dr 


dr 


dz 


}_(dTi  \_dn  PTl 

Pr  1 dr 2 ^ r dr  dz 2 


and  in  the  solid 


dTs  = 1 fdTs  1 dTs  d2Ts 

dr  Ps  1 dr2  r dr  dz 2 


(9) 


where  Pr  = v / kl  is  the  Prandtl  number  of  the  melt,  and  Ps  = v / k s',  here  kl,  k,s  are  the 
thermal  diffusivities  of  the  liquid  and  solid  phases  respectively. 

The  solid  temperature  is  prescribed  on  the  inner  rigid  wall,  r = r o,  and  on  the  outer 
rigid  wall,  r = r2,  we  have 


u = w = 0,  Tl  = 1,  v = (I/77)  COST. 
On  the  interface  r = 77 / ( 1 — 77)  -f  h,  we  have 

u = w = 0, 

v = \ 1 H —h)  cos r, 

V V J 


Tl  = Ts  = 7/C, 


(10) 

(11) 

(12) 

(13) 


and 


dh  _dTL_dhdTL_  ( 9TS  dhdTs\ 

dr  dr  dz  dz  ^ ^ dr  dz  dz  ) 


(14) 


Here  K is  twice  the  mean  curvature  in  cylindrical  coordinates,  7 = TmF/(LAT)  is 
proportional  to  the  surface  energy,  q = ks/ki,  is  the  ratio  of  the  thermal  conductivities  in 
the  crystal  and  the  melt  respectively,  and  C — uLv/(ki,AT ) is  proportional  to  the  latent 
heat  Lv  released  upon  solidification. 

The  Taylor  number  may  be  related  to  the  Reynolds  number  Re  used  in  previous  work 
[25]  by 

9n 

(15) 


Ta  = -r— — Re2, 


1 -77 
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where  Re  = SIL2 /v.  This  definition  of  the  Taylor  number  is  identical  to  that  used  by 
Hall  et  al  [26,  27]. 


3.1.1.  Base  State 


We  seek  a base  state  that  is  only  dependent  on  r and  t.  The  base  state  flow  field  is  given 
by  — u/°)  = 0 and  = Real[Vo(^)eIT],  and  where  Vo(r)  satisfies 


dV0  ldV0 

dr2  ^ r dr 


- (p  + iu)  v°  - °> 


(16) 


subject  to  Vo(7T)  = 1 and  Vo(r2)  = 1/rj.  Though  the  solution  may  be  written  in  terms  of 
Kelvin  functions,  we  choose  to  integrate  it  numerically  at  the  collocation  points  used  in 
the  subsequent  discretization.  The  thermal  fields  satisfy  ordinary  differential  equations 
in  r;  the  time-independent  base  states  for  the  temperature  fields  are  given  by 


and 


ln[r(l  -7})/ri\ 
In  7] 


(17) 


In  Hi  -’?)/’>] 

gin  77 

The  base  state  interface  shape  is  given  by  — 0.  The  base  state  pressure  gradient  may 
be  determined  from  the  radial  momentum  equation. 
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3.1.2.  Linearized  Equations 


We  perturb  the  base  state  and  denote  the  disturbance  quantities  by  a superscript  (1). 
Assuming  normal-mode,  small-amplitude  disturbances,  we  write 


0 1 

V 

?/0)(r,  r) 

u(1)(r,  r) 

w 

0 

w^\r,  r) 

p 

— 

PIO,(r,T) 

+ 

p(1)(r,  r ) 

Tl 

1 fV) 

^L)(r.r) 

Ts 

1 f\r) 

^1}(r,T) 

V h ) 

i 0 / 

1 k(1)M  j 

exp(iaz), 


(19) 


where  a is  the  axial  wavenumber.  Eliminating  the  pressure  and  axial  velocity  w^\ 
the  linearized  equations  become 

a,  {DD.  - a2)  ^ + a2^Ta^V)  - {D2  D]  - 2a2 DD.  + a4)u^\  (20) 


dvM 


a f = (DD*  — a2) 


dr 

on  ri  < r < r2,  and  in  the  solid  r0  < r < r\, 

dT{s] 

dr 


+ Z)(7l°))uW  = ±-(D,D-a1)Tf\ 


(21) 

(22) 


UJ- 


= rs{D’D-“2)Tp- 


(23) 


Here  D = d/dr  and  Z)*  = d/dr  -f  1/r. 

At  the  inner  solid  wall  r0  = t?s/(1  — v),  we  set  — 0,  and  at  the  outer  wall 
r = r2  = 1/(1  — 77),  we  require 


= DuW  = T™  = 0. 

At  r —r  1 = 77 / ( 1 — T])  the  crystal- melt  interface  boundary  conditions  are 

v,M  = DuW  = 0, 


(24) 


(25) 
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v™  + 


i - V 
V 


cos  r 


+ DvW 


h(1)  = 0, 


+ DT^  - qDT^  = 0, 

+ (Dlf>)hi''>  = Tlsl)  + (Blf>)h(l'>/q  = 7 (l/rj  - a2)  hll>. 


(26) 

(27) 

(28) 


The  above  set  of  equations  have  temporally-periodic  coefficients  that  vary  spatially 
as  well;  the  solutions  may  be  determined  from  an  application  of  Floquet  theory. 


3.2.  Narrow  Gap 


3.2.1.  Scalings  and  Base  State 


We  shall  also  study  the  simplified  system  of  equations  that  results  in  the  narrow  gap 
limit.  All  scales  remain  the  same  except  that  we  now  scale  the  radial  coordinate  as 
r'  — Ri  + Lx ; in  nondimensional  form,  r = r\  + x.  The  melt/solid  interface  that  is  at 
ri  =77/(1  — 77)  in  the  base  state  will  be  located  at 

r = ri  + h,  (29) 

where  the  nondimensional  interface  amplitude  h is  scaled  with  the  gap  width.  In  nondi- 
mensional form,  the  domain  becomes  Xo  < x < 1,  where 

£0  = (77s  — 77)/(l  — 77)  (30) 


with  the  melt  occupying  h < x < 1.  In  the  linearized  theory,  the  melt /solid  interface 
location  h will  be  taken  close  to  zero. 

In  this  limit,  we  assume  1 — rj  = e <C  1,  and  that  the  Taylor  number  and  the  dependent 
variables  remain  0(1).  We  then  have  = x + 0(e),  T ^ = x/q  + O(e),  and  the  base 
state  velocity  profile  now  satisfies  u — u/°)  = 0 and 


v(°\0 , r)  = i/°)(1,t)  = cost 

The  solution  for  the  azimuthal  component  is 

. cosh  [\/iw(i  — l/2)j  1 


= Real  { e*T- 


cosh(%/iu>/2) 


(31) 


(32) 
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The  pressure  p may  be  again  determined  from  the  radial  momentum  equation. 

3.2.2.  Linearized  Equations 

We  may  disturb  the  basic  state  as  in  Eq.  (19),  and  eliminating  the  pressure  and  axial 

\ 

velocity  w^\  the  linearized  equations  become 

( d2 


LJ 


dx2 


— a 


2 ^ ^ U ^ ^ i 2m  (0)  (1) 

1 +o  larV^ 


dr 


u- 


dv^  dv^ 
dr  + dx 

dT{rl) 


V1)  = 


UJ- 


dr 


on  0 < x < 1,  and  in  the  solid  x0  < x < 0, 


d2 

dx2 

' d2 

dx2 


— a \ uK 


- a2 1 «<». 


+ «W  = 


1 (£.  _ aA  T« 

Pr  l dx2  L ’ 


(33) 

(34) 

(35) 


UJ- 


dT{s1] 

dr 


At  the  inner  solid  wall  x0  = (775  — 77 )/( 1 — 77), 

T M = 0. 

At  the  outer  wall  x = 1, 


_L  (£.  _ T« 

Ps  Uz2  s ' 


(36) 


(37) 


v<»  = „(')  = ^ = I'M  = 0. 

dx 


At  the  crystal-melt  interface  referred  to  x = 0, 


(38) 


ID  du(1)  „ 

U(1)  = -T—  = 0, 
ox 

(39) 

+ afV>  = o, 

dx 

(40) 

r dhM  dT^]  dT{s1] 

U dr  dx  q dx  ’ 

(41) 

+ A(1)  = + -h™  = -ia2hW. 

q 

(42) 

Because  of  the  spatial  and  temporal  variation  in  we  must  solve  this  problem 

numerically  for  most  values  of  the  parameters. 
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4.  Numerical  Formulations 

In  general,  the  above  set  of  equations  must  be  solved  numerically.  We  have  carried  out 
a pseudospectral  discretization  in  the  spatial  variable  r;  the  resulting  system  of  ordinary 
differential  equations  subject  to  algebraic  constraints  may  be  integrated  over  one  period 
in  time,  and  the  Floquet  exponents  computed;  we  may  then  determine  the  stability  of  the 
base  state.  We  may  also  Fourier  expand  the  solutions  in  time  and  find  the  growth  rate 
of  each  mode  by  solving  the  resulting  two-point  boundary  value  problems  as  described 
elsewhere [25].  We  have  employed  both  procedures  as  a check  on  our  numerics. 


4.1.  Wide  Gap 


We  discretize  the  derivatives  in  the  radial  direction  using  the  standard  Chebyshev  pseu- 
dospectral method  (see,  e.g.,  [37,  38]).  It  is  convenient  for  this  method  to  rescale  the 
radial  variable.  In  the  liquid  we  set  r = r i + |(£  + 1)  for  — 1 < £ < 1,  and  in  the  solid 
we  set  r = r0  -f  |(£  + l)x0.  Then  in  the  liquid  we  have  d/dr  = 2 d/d£  and  in  the  solid 
we  have  d/dr  = —(2 /x0)d/d£i. 

We  use  the  points  = cos  jir /n  for  j = 0, 1, . . . , n,  so  that  the  interface  is  located  at 
£o  = 1-  At  the  collocation  points  we  use  the  Chebyshev  derivative  matrix  D [38],  which 
has  the  property  that  at  the  collocation  points  (j  the  derivative  g'-  of  an  n-th  degree 
polynomial  g{ £)  is  given  exactly  in  terms  of  its  collocation  values  g k by  the  expression 


DjkQk- 


k—0 


(43) 


Higher  derivatives  are  represented  by  powers  Dm  of  the  matrix  D.  We  write  g'  = Dj^gk 
and  thereby  let  the  sum  over  the  repeated  index  k be  implied.  Dij  is  the  differentiation 
matrix  in  the  melt,  and  Dtj  is  the  differentiation  matrix  in  the  solid. 

For  convenience,  we  drop  the  superscript  on  the  disturbance  quantities;  at  the  interior 
points  in  the  melt  (j,  j — 1, . . . , n — 1,  we  have  the  discretized  equations 


^§7  + ^Dtj  + ^7r*K'0)] 


Dij  + Dij/ n - Sij(a 2 + l/rt2) 


J31 


(44) 


-14- 


MODULATED  TAYLOR-COUETTE  FLOWS 


November  4,  1992 


a;[^  + AJ/r1-^(a2  + l/rt2)j-^  + 


duj  a2r] 


,(o) 


(i  _ „)Ta7-Mj  = W + 2D'’,r'~ 


LO- 


(3/r2  + 2 a2)D2tj  + (3 /rf  - 2a2/r,)D„  - (3/r4  - 2a2/r2  - a4)S,j 

dTu  dT (!> 


u 


31 


dr  dr 
and  in  the  sohd,  we  have 


r)rTvJ)  1 

+ = — {D%  + Dij/ti  - Saa1)  TLj] 


(45) 

(46) 


"fr  = i (4  + 6.,/r.  - «,,a2)  TS, 


Here  8tJ  is  the  Kronecker  delta,  and  we  sum  over  the  repeated  index  j. 

In  both  phases,  the  subscript  0 indicates  the  rigid  boundary  and  the  subscript  n 
indicates  the  interface.  We  also  have  a differential  equation  for  the  interface  amplitude 
given  by 

dh 

~ ~DnjTLj  + qDnjTsj.  (47) 

Instead  of  satisfying  the  differential  equation  for  u at  X\  and  xn_i,  we  apply  the  algebraic 
conditions 


DijUj  — 0 and  Dnjiij  = 0. 


(48) 


The  boundary  condition  on  the  azimuthal  velocity  becomes  the  algebraic  equation 

^ DnjVj — — cos t ^ = 0.  (49) 


From  continuity  of  the  temperature  field  at  the  interface  we  have  the  final  two  algebraic 
conditions 


and 


1/r2  J h ^ 


0, 


(50) 


1/r2  J 


0. 


(51) 


There  are  then  An  — 5 differential  equations  and  5 algebraic  equations  to  determine  the 
An  unknowns.  The  computer  code  DASSL  [36]  is  used  to  solve  the  differential-algebraic 


-15- 


MODULATED  TAYLOR-COUETTE  FLOWS 


November  4,  1992 


system  in  time,  in  a manner  similar  to  that  described  by  Murray  et  al. [34].  The  Floquet 
analysis  is  performed  by  writing  the  solution  vector  F in  the  form 

F{r,t)  = f(r,t)eat, 

where  the  amplitude  /(r,  t ) is  periodic  in  time,  and  the  complex  growth  rate  a determines 
the  linear  stability  of  the  system.  The  numerical  determination  of  a is  implemented  by 
constructing  the  K x K fundamental  solution  matrix,  where  K = 4n  — 5 is  the  number 
of  differential  equations.  The  eigenvalues  of  this  matrix  are  the  Floquet  multipliers  from 
which  a is  obtained. 

4.2.  Narrow  Gap 

We  also  solve  the  narrow  gap  problem  by  using  a pseudospectral  discretization  in  the 
spatial  variables,  and  this  reduces  the  problem  to  ordinary  differential  equations  in  time 
coupled  to  algebraic  conditions  at  the  boundary.  The  discretized  equations  for  the  ve- 
locity components  become 


" K - °2^)  fr  = 

(D%  - 2a2 Dl  + a45„)  u,  - a'Tavj0’^, 

(52) 

dvl 

UJ = 

dr 

(Dij  ~ sija2)  vj  - DikV^SijUj, 

(53) 

The  equation  for  the  temperature 

in  the  melt  becomes 

8TQ 

“ dr 

s 

[d2,  - a2StJ ) tW  - u„ 

(54) 

and  in  the  solid, 

dTf> 

U dr 

1 

Ps 

— D2  — a2ti 

a otJ  iSj. 

x0 

(55) 

In  both  phases,  the  subscript  0 indicates  the  rigid  boundary  and  the  subscript  n indicates 
the  interface.  The  final  differential  equation  determines  the  interface  location 


+ Dn,T$  - 9i-0n,T$  = 0.  (56) 
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Instead  of  satisfying  the  differential  equation  for  u at  Xi  and  xn_i,  we  apply  the  algebraic 
conditions 

Di juj  = 0 and  DnjUj  = 0.  (57) 

The  boundary  condition  on  the  azimuthal  velocity  becomes  the  algebraic  equation 

dW  + DnjvfhW  = 0.  (58) 

From  continuity  of  the  temperature  field  at  the  interface  we  have  the  final  two  algebraic 
equations 

T&>  + /i(1)  = t£>  + q-'h^  = -7a2fc<l>.  (59) 

We  again  have  4n  — 5 coupled  differential  equations  and  5 algebraic  conditions  to  deter- 
mine the  4n  unknowns. 


5.  Numerical  Results 

First,  convergence  of  the  results  for  the  wide  gap  model  to  those  for  the  narrow  gap 


Wide 

Gap 

V 

Ta 

0.69 

991.23 

0.9469 

901.26 

0.9832 

899.02 

0.9917 

897.06 

Table  1:  Linear  theory  results  verifying  convergence  of  wide  and  narrow  gap  codes 
for  uj  = 28.9  and  a = 1.866.  Note  that  the  result  in  the  first  row  is  in  very  good 
agreement  with  previous  results [25];  there,  Ta  = 991.21  for  the  same  parameters 
except  7 = 0.  The  narrow  gap  result  (rj  — 1)  is  Ta  = 895.07.  Here  the  crystal 
width  is  fixed  at  rx  — r0  = 1.33  and  rjs  — 7 — (77  — r0)(l  — 77). 


model  can  be  demonstrated,  as  is  illustrated  in  Table  1. 

Neutral  curves  for  several  values  of  the  Prandtl  number  for  the  wide  gap  model  are 
shown  in  Figure  3.  For  small  Prandtl  numbers,  the  neutral  curve  approaches  that  of 
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the  rigid-walled  case,  while  for  increasing  Pr,  we  see  a destabilization  of  the  purely 
hydrodynamic  mode  and  an  accompanying  decrease  in  the  critical  wave  number.  The 
critical  Taylor  number  becomes  approximately  proportional  to  Pr-1  as  is  shown  in  Figure 
4 and  Table  2.  As  the  Prandtl  number  increases  through  unity,  the  neutral  curve  rapidly 
approaches  the  limiting  shape  denoted  by  Pr  = oo.  This  limiting  neutral  curve  is  the 
result  of  an  asymptotic  analysis  for  large  Prandtl  number  which  is  discussed  in  the 
following  section. 

The  destabilization  of  the  centrifugal  instability  occurs  precipitously  as  the  Prandtl 
number  is  increased.  This  destabilization  is  apparently  caused  by  the  strong  coupling  of 
the  temperature  field  with  the  flow  in  the  melt  and  the  subsequent  transport  of  heat  from 
the  interface.  When  Pr  1,  the  dominant  mode  of  heat  transfer  is  from  convection, 
and  the  temperature  field  is  strongly  affected  by  the  disturbance  flow.  For  Pr  *C  1,  heat 
transfer  is  dominated  by  conduction,  and  the  flow  disturbance  discussed  here  has  little 
effect  on  the  interface  shape [22]. 

Examination  of  the  Fourier  coefficients  of  the  unstable  modes  shows  that  the  response 
is  such  that  the  azimuthal  velocity  consists  of  only  odd  temporal  harmonics  and  the 
remaining  dependent  variables  consist  of  only  even  temporal  harmonics;  this  type  of 
response  was  termed  “Type  I”  in  [25],  where  modes  having  a different  symmetry  were 
also  obtained  in  the  rigid-walled  system.  It  may  further  be  seen  that  the  solutions  are 
dominated  by  the  mean  (time  independent)  and  fundamental  components  (the  lowest 
frequency  terms  from  the  Type  I solutions). 

6.  Asymptotic  theory  for  Pr  1 

We  next  describe  an  asymptotic  analysis  in  the  limit  Pr  — ■>  oo  for  the  destabilization  of 
the  centrifugal  instability  in  the  two-phase  system.  We  describe  the  narrow  gap  results 
in  detail  for  ease  of  discussion,  and  simply  report  the  results  for  the  wide  gap  case. 
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6.1.  Narrow  Gap 

If  the  Prandtl  number  is  large  enough,  limiting  behavior  occurs  in  the  solutions  to  the 
linear  problem.  Some  numerical  results  for  several  values  of  the  Prandtl  number  with 
all  other  parameter  values  fixed  are  given  in  Table  2.  The  other  parameters  are  fixed 
at  values  appropriate  for  the  experiments  with  SCN.  We  note  that  the  simpler  problem 
with  steady  rotation  was  also  amenable  to  large  Prandtl  number  asymptotics  [22].  We 


Pr  = Ps 

Ta 

T =PrTa 

1.0 

17315. 

17315. 

22.8 

899.79 

20515. 

100. 

206.44 

20644. 

1000. 

20.678 

20678. 

Table  2:  Results  of  Floquet  theory  with  pseudospectral  spatial  discretization  for 
the  narrow  gap  limit  for  increasing  Pr  number.  For  all  of  the  values  in  the  table, 
uj  = 28.9,  a = 1.75,7  — 6.2  x 10~6,  C = 409,  and  q = 1.009.  Note  that  the 
product  PrTa  tends  to  a constant  as  the  Prandtl  number  increases. 


now  carry  out  an  asymptotic  analysis  for  the  modulated  problem  in  the  large  Prandtl 
number  limit. 

One  approach  to  use  at  this  point  is  Floquet  theory  with  Fourier  analysis,  as  discussed 
in  [25].  This  is  suggested  by  an  examination  of  the  equations,  with  due  regard  to  the 
occurrence  of  harmonics.  As  in  Hall  et  al. [26,  27],  one  expects  the  Fourier  coefficients  to 
depend  on  Pr.  The  Fourier  decomposition  approach  shows,  and  it  is  confirmed  by  exam- 
ination of  the  numerically-determined  Fourier  coefficients  of  the  pseudospectral  solution, 
that  the  following  expansions  should  be  used: 


■j/1)  ^ 


u1(x)etT  + ^^e3lT  + ...  + c.c. 


Pr 


1 


u(1)  ~ ju0(x)  + 


H*V'r  + -4,T 


e™T  + . . . + c.c. 


if1’  ~ f0(x)  + 


Pr 


(60) 

(61) 

(62) 
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if  > ~ fso(x)  + 

h ~ h0  + 


?S2(*)  2*t 


Pr 


e + . . . + c.c. 


^-e2lT  + . . . + c.c. 
Pr 


(63) 

(64) 


and  that  the  Taylor  number  should  be  scaled  as  Ta  = T /Pr.  Here  we  have  assumed 
scalings  for  the  first  harmonics  of  the  temperature  in  the  solid  and  the  interface  suggested 
by  the  interfacial  temperature  condition;  we  defer  to  the  appendix  further  discussion  of 
the  scaling  behavior. 

The  problem  for  the  mean  and  fundamental  coefficients  becomes 


iouv  i = ( D 2 — a2)v  i, 

(65) 

(. D2  — a2)2u0=  ±a2T{V0vl  + c.c.), 

(66) 

uq  = ( D 2 — a2)T0, 

(67) 

on  0 < x < 1 and  on  x0  < x < 1 we  have 


(. D 2 - a2)TS0  = 0.  (68) 

At  the  rigid  boundaries  all  of  the  dependent  variables  are  zero,  and  at  the  interface  we 
must  have 


Uq  - Duq  = v1  + \DV0h0  — 0, 

(69) 

To  + ho  = Tso  + Q % — 

(70) 

DTq  — qDTso  — 0. 

(71) 

Here  D = d/dx.  This  time-independent  system  was  solved  for  numerically  using  SU- 
PORT  [39]  from  the  SLATEC  Common  Math  Library  [40];  this  leading-order  boundary 
value  problem  represents  a drastic  simplification  from  the  numerical  implementations  of 
Floquet  theory  described  earlier  in  this  paper.  The  problem  may  also  be  solved  analyt- 
ically to  obtain  a complete  closed  form  solution.  However,  this  closed  form  solution  is 
rather  complicated  and  no  advantage  is  gained  in  interpreting  results.  The  solution  can 


-20- 


MODULATED  TAYLOR-COUETTE  FLOWS 


November  4,  1992 


be  expressed  as  the  sum  of  a complementary  part,  which  is  independent  of  the  Taylor 
number,  and  a particular  solution  which  is  proportional  to  the  Taylor  number. 

We  choose  to  solve  the  complementary  part  of  the  problem  in  closed  form,  and  then 
determine  numerically  the  remaining  particular  part  to  determine  T — PrTa.  We  find 

--DVo(O)  sinh  \fa?  + zu;(l  — x ) 

2 sinh  \/a2  + iu> 

u0(x)  = Tu0(x),  (73) 

- . . . sinh  a(l  — x)  _ ~ . 

r„(a0  = -(l+7<i2) At >-  + TT0{x),  (74) 

smn  a 

and 


Ui(x)  = 


TSo{x ) = - Q + 7a2 


sinh  a(xo  — x) 


sinh  axo 

For  the  particular  solutions  (denoted  by  a ~),  we  then  have 


( D 2 - a2)2u0  = \a2(VoV*x  + c.c.), 
u0  = ( D 2 - a2)T0 , 

subject  to 


(75) 


(76) 

(77) 


uo(0)  = Du o(0)  = 7o(0)  = tio(l)  = Biio(l)  = T0(l)  = 0.  (78) 

Once  these  solutions  are  known,  we  may  then  determine  the  Taylor  number  via  Eq.  (71). 

In  the  next  section  we  present  the  corresponding  simplified  problem  for  the  wide  gap 
case,  and  then  provide  a discussion  of  the  asymptotic  results  for  both  the  narrow  and 
wide  gap  cases. 

6.1.1.  Wide  Gap 

We  employ  the  same  scales  as  in  the  narrow  gap  case  to  obtain  equations  for  the  wide  gap 
case  that  are  valid  for  large  Pr.  The  problem  for  the  mean  and  fundamental  coefficients 
becomes 

iujv  i = ( DD * — a2)vi,  (79) 


-21- 


MODULATED  TAYLOR-COUETTE  FLOWS 


November  4,  1992 


(D2D2  - 2 a2DD,  + a4)i0  = ^T-(V0v\  + c.c.), 

2 r 

(°) 

— = (D.D  - a2)%, 
or 


on  0 < x < 1,  and  on  Xq  < x < 1 we  have 


(80) 

(81) 


( D.D  - a2)Tso  = 0.  (82) 

At  the  rigid  boundaries  all  of  the  variables  are  zero,  and  at  the  interface  we  must  have 
u0  = Duq  = Vi  + \ ( DV0  - l/n)  h0  = 0,  (83) 

qji(o)  _ _ dT^  - 

To  ~\ rr—ho  = Tso~\  rr—  ho/q=  l/r\  — a2  'yho,  (84) 

or  Or  1 * 

DTo  - gDTso  = 0.  (85) 


Note  that  dT^/dr  = — l/(rln77).  This  system  is  again  solved  as  a boundary  value 
problem  with  SUPORT;  the  asymptotics  have  eliminated  the  need  for  Floquet  theory. 
From  the  standpoint  of  the  design  of  the  experiment,  this  is  a big  advantage,  as  we  shall 
see. 

Part  of  the  wide  gap  problem  may  be  solved  in  closed  form  in  terms  of  modified 
Bessel  functions;  however,  for  convenience  we  choose  to  find  the  solutions  numerically 
using  SUPORT.  We  discuss  the  results  of  the  asymptotic  analysis  in  the  following  section. 

7.  Results  and  Discussion 

Numerical  results  from  the  asymptotic  approximations  Eq.  (65)-(68)  and  Eq.  (69)-(71) 
for  the  narrow  gap  equations  and  from  Eq.  (79)-(82)  and  Eq.  (83)-(85)  for  the  wide  gap 
equations  are  displayed  as  curves  labeled  “Pr  = oo”  in  Figures  4 and  5.  These  results  are 
for  a fixed  frequency  of  uj  = 28.9.  The  neutral  curves  from  the  full  Floquet  theory  rapidly 
approach  these  limiting  results  for  Pr  > U particularly  for  smaller  wavenumbers  in  the 
neutral  curve.  The  entire  neutral  curve  from  the  asymptotic  analysis  can  be  obtained  in 
two  orders  of  magnitude  less  processing  time  than  that  required  to  compute  a single  data 
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point  by  using  the  full  Floquet  theory  implementation.  The  limiting  behavior  where 
critical  wavenumber  is  decreased  and  the  Taylor  number  scales  like  the  inverse  of  the 
Prandtl  number[25]  occurs  even  as  the  Prandtl  number  is  near  unity,  and  it  is  a good 
approximation  to  the  experimental  case  of  Pr  = 22.8  for  SCN. 

Figure  5 displays  the  limiting  behavior  for  the  narrow  gap  approximation.  The  smaller 
wavenumber  side  of  the  curve  is  virtually  unaffected  by  the  Prandtl  number  variation,  and 
as  the  Prandtl  number  increases,  the  higher  wavenumber  side  approaches  the  asymptotic 
limit  (Pr  = oo). 

The  wide  gap  and  narrow  gap  results  are  compared  in  Figure  6;  the  results  become 
insensitive  to  the  gap  width  once  77  is  larger  than  about  0.69.  As  was  found  in  the 
steadily-rotating  case  [23],  the  narrow  gap  results  bound  the  wide  gap  results  from  below. 
It  appears  reasonable  then  to  use  the  narrow  gap  results  to  design  experiments  in  the 
large  Pr  regime. 

Table  3 displays  asymptotic  and  numerical  results  for  the  wide  gap  model;  in  this 
table,  the  critical  wave  number  is  found  as  a function  of  the  frequency  for  each  method. 
It  is  clear  from  the  table  that  as  long  as  the  frequency  is  not  too  low,  the  asymptotic 
results  do  a good  job  of  approximating  the  Floquet  theory  results.  The  numerical  results 
are  in  qualitative  agreement  with  the  preliminary  experimental  results.  For  example, 
the  wavenumber  of  the  instability  from  the  experiment  is  about  1.7  for  the  conditions 
uj  = 28.9,  rj  = 0.69,  and  775  = 0.286.  This  is  in  good  agreement  with  the  theoretically 
determined  value  of  a — 1.87.  The  critical  Taylor  number  from  the  asymptotic  theory 
is  Ta  = 987,  while  the  system  is  observed  to  be  unstable  in  the  preliminary  experiments 
for  this  frequency  at  Ta  = 2038.  This  instability  is  still  much  lower  than  the  marginal 
theoretical  value  of  Ta  = 93857  at  a critical  wavenumber  of  a = 5.059  for  the  rigid- 
walled  system  in  which  there  is  no  crystal- melt  interface[25];  there  is  no  experimental 
data  available  for  the  rigid-walled  system  to  our  knowledge. 

The  asymptotic  results  above  have  been  found  for  a nondimensional  frequency  that 
is  0(  1)  with  respect  to  the  Prandtl  number.  We  find,  however,  that  agreement  with 
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Asymptotic 

Floquet 

u> 

ac 

PrTa 

ac 

PrTa 

5.0 

1.449 

23193.9 

1.241 

30442.0 

7.07 

1.491 

15868.9 

1.387 

17993.5 

10.0 

1.557 

13075.9 

1.502 

13880.8 

15.0 

1.662 

13616.7 

1.637 

13961.4 

21.0 

1.765 

16915.9 

1.757 

17110.5 

28.9 

1.868 

22511.4 

1.866 

22599.6 

40.0 

1.969 

30032.3 

1.978 

30001.0 

50.0 

2.032 

35529.0 

2.043 

35425.6 

60.0 

2.078 

39996.4 

2.086 

39854.8 

70.0 

2.111 

43846.2 

2.120 

43687.9 

80.0 

2.137 

47385.1 

2.145 

47220.2 

90.0 

2.158 

50790.7 

2.166 

50624.5 

100.0 

2.175 

54152.8 

2.182 

53985.7 

110.0 

2.190 

57510.5 

2.197 

57343.7 

Table  3:  Results  of  Floquet  theory  and  asymptotic  theory  for  the  wide  gap  case. 
In  this  case,  the  critical  wave  number  is  found  to  the  nearest  10-3  for  a fixed 
frequency.  The  parameters  used  for  this  table  are  C = q = 1.0,  Pr  = Ps  = 
22.8,  7 = 0,  77  = 0.69,  775  = 0.286. 
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high  frequency  results  from  the  full  implementation  of  Floquet  theory  is  remarkably 
good.  Table  3 and  Figure  7 illustrate  these  results.  Although  the  high  frequency  results 
agree  rather  well,  the  agreement  in  the  low  frequency  cases  is  not  so  good.  We  believe 
that  this  occurs  because  in  the  high  frequency  limit  the  time-independent  modes  of  the 
solution  again  dominate[34],  in  a manner  similar  to  that  just  demonstrated  for  the  high 
Prandtl  number  regime.  Though  we  have  not  carried  out  such  an  analysis,  the  results 
are  reminiscent  of  those  that  would  occur  in  an  averaging  analysis  at  high  frequency. 
In  the  low  frequency  regime,  on  the  other  hand,  the  temporal  behavior  is  complicated 
and  cannot  be  described  by  a single  mode  in  time[34].  The  computation  at  the  lowest 
frequencies  is  difficult  for  the  pseudospectral  approach  with  Floquet  theory,  and  the 
computations  are  best  performed  there  with  the  Fourier  approach  of  [25]. 

Although  the  Floquet  theory  and  asymptotic  approximations  agree  closely,  prelim- 
inary experimental  results  are  only  in  qualitative  agreement  with  either  theory.  The 
instability  has  been  observed  only  above  the  neutral  curves  given  by  either  theory  (see 
Figure  4 in  [25]).  The  effect  of  the  deformable  crystal-melt  interface  is  pronounced;  both 
the  theoretical  and  experimental  results  show  that  the  two-phase  system  is  significantly 
less  stable  than  the  rigid- walled  system. 

8.  Conclusion 

We  have  modeled  the  destabilization  of  the  purely  azimuthal  base  flow  in  the  Taylor- 
Couette  geometry  with  a crystalline  inner  cylinder  under  pure  torsional  modulation. 
We  have  developed  a pseudospectral  discretization  in  the  spatial  variables  and  found 
the  eigenvalues  of  the  fundamental  solution  matrix  of  the  resulting  differential/algebraic 
system  in  order  to  determine  stability  or  instability [34]. 

In  addition,  we  have  been  able  to  carry  out  asymptotic  analyses  in  the  limiting  case 
of  large  Prandtl  number  in  the  Taylor-Couette  geometry  with  a crystalline  inner  cylinder 
under  pure  torsional  modulation.  In  an  analysis  similar  to  the  work  of  Hall  et  aZ. [26 , 27] 
we  have  carried  out  a Fourier  analysis  with  the  amplitude  of  each  mode  depending  on 
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the  large  parameter,  in  this  case  the  Prandtl  number.  Because  the  leading-order  mode 
is  primarily  time-independent,  the  “full”  numerical  analysis  of  the  problem  involving 
Floquet  theory  is  simplified  to  a single  boundary  value  problem  at  leading  order.  This 
results  in  a great  reduction  in  computational  effort.  Because  the  coefficients  of  the  higher- 
order  modes  also  display  boundary-layer  structure  in  space  (see  appendix),  they  appear 
to  be  even  smaller  than  the  expected  scalings;  this  helps  make  the  analysis  a particularly 
effective  approximation.  The  approximation  remains  good  in  the  high-frequency  regime 
where  a time-independent  (or  steady-streaming)  mode  of  the  solution  components  is 
dominant. 

Given  the  preliminary  nature  of  the  experiments  performed  with  succinonitrile,  the 
agreement  in  the  onset  of  the  instability  seems  good.  The  wavenumber  is  reasonably 
approximated  by  the  linear  theory  value.  It  appears  that  in  the  large  Prandtl  number 
regime,  the  transport  of  heat  is  strongly  affected  by  the  disturbance  flow  of  the  melt 
and  that  the  flow  in  the  melt  is  destabilized  by  the  deformable  interface.  More  extensive 
experiments  would  serve  to  clarify  the  correspondence  between  the  experimental  and 
theoretical  results. 

The  general  question  of  predicting  whether  a given  hydrodynamic  instability  will  be 
strongly  modified  if  a rigid  boundary  is  replaced  by  a deformable  crystal-melt  interface 
is  still  poorly  understood.  Several  studies  have  shown  that  instabilities  associated  with 
critical  layers  in  the  interior  of  the  fluid  are  not  affected  strongly  by  the  presence  of  the 
crystal-melt  interface.  This  is  not  unexpected,  since  for  these  instabilities  the  perturba- 
tions to  the  flow  are  generally  confined  to  the  immediate  vicinity  of  the  critical  layer,  and 
influence  of  the  boundary  conditions  of  the  perturbed  flow  is  not  large.  Examples  of  this 
type  include  the  instabilities  associated  with  Poiseille  flow  [20],  the  asymptotic  suction 
profile  [16],  and  the  buoyant  instability  of  the  parallel  flow  inside  a vertical  annulus  with 
lateral  heating  [21,  20]. 

For  instabilities  of  a less  local  nature,  for  which  the  relevant  length  scale  of  the  instabil- 
ity is  the  container  width,  the  possibility  of  flow/interface  interaction  is  correspondingly 
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greater.  Examples  of  this  type  include  the  shear  instability  of  the  parallel  flow  inside  a 
vertical  annulus  with  lateral  heating  [21,  20],  Rayleigh-Benard  convection  [17],  steady 
Taylor- Couette  flow  [22,  23],  and  the  modulated  Taylor- Couette  flow  considered  here. 
For  these  flows,  strong  couplings  are  observed  for  large  Prandtl  numbers,  except  in  the 
case  of  Rayleigh-Benard  convection,  where  the  modification  of  the  linear  stability  of  the 
system  is  not  found  to  be  large  [17].  A common  feature  of  the  instabilities  which  are 
strongly  affected  by  the  crystal-melt  interface  is  the  presence  of  shear  in  the  base  state. 
With  shear,  the  linearized  no-slip  boundary  condition  for  the  perturbed  velocity  pro- 
duces a direct  coupling  between  the  tangential  perturbation  velocity  and  the  interface 
deformation;  this  coupling  is  absent  for  a quiescent  base  state,  as  in  Rayleigh-Benard 
convection. 
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Appendix  A.  Higher  Order  Modes 

We  now  discuss  why  the  lowest-order  terms  dominate  for  Pr  = Ps>>lby  examining  the 
problem  for  the  next-order  nonzero  harmonics.  As  mentioned  previously,  the  scalings  for 
the  harmonics  is  assumed  from  the  form  of  the  boundary  conditions.  From  the  numerical 
results,  however,  it  can  be  seen  that  the  scales  for  the  coefficients  of  the  harmonics  for 
the  interface  position  and  the  solid  temperature  field  are  of  higher  order  than  expected. 
The  problem  for  the  next-order  harmonic  Fourier  coefficients  is 

3ivv3  = (D2  — a2)v3  - ±(DV0u2  + ^DVqU4),  (Al) 
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2 iw(D2  - a2) u2  = ( D 2 - a2)252  - \a2T{V0v^  + (A2) 

2iu>T2  + U2  = — (-D2  — a2)-^2)  (A3) 

on  0 < x < 1 and  on  xq  < x < 1 we  have 


2 iwfS2  - A(D2  - a2)fS2  = o.  (A4) 

Ps 

At  the  rigid  boundaries  similar  homogeneous  conditions  hold,  and  at  the  interface  we 
have 


u2  — Du2  — u3  + 2DV0h2  + 2pr^^o  ^4  — 0, 

(A5) 

T2  + h2  = Ts2  + gsh  2 — — 7fl2^2) 

(A6) 

2iu) C,h2  T DT2  — gDTs2  — 0. 

(A7) 

Here  the  asterisk  denotes  complex  conjugation.  In  order  to  write  these  equations  we  have 
assumed  that  the  first  harmonic  coefficient  for  the  interface  shape  and  the  temperature 
in  the  solid  scale  hke  1/Pr;  this  is  suggested  by  the  interfacial  temperature  conditions. 
Now  to  leading  order  for  Pr  1 in  these  Fourier  coefficients  some  of  the  terms  in  the 
equations  drop  out.  In  particular,  the  heat  equation  in  the  melt  yields 

T2  — -zr—u2]  (A8) 

zzu; 

since  u2  = Du2  — 0 at  both  the  rigid  wall  and  the  interface,  then  T2  = DT2  — 0 there  as 
well.  There  will  be  an  0(1/Pr)  correction  to  T2.  The  heat  equation  in  the  solid  yields 
that  Ts2  — 0.  The  interfacial  temperature  conditions  then  imply  that  h2  = 0 to  leading 
order.  Because  DT2  = 0,  the  heat  flux  at  the  interface  vanishes  to  leading  order  as  well. 
The  Fourier  coefficients  for  the  leading  order  temperature  in  the  solid  Ts2  appears  to 
be  zero  by  a similar  process.  The  velocity  components  u2  and  u3  must  come  from  the 
numerical  solution  of  the  coupled  boundary  value  problem  given  by  Eq.  (A1)-(A4)  with 
boundary  conditions  Eq.  (A5)-(A7). 

The  problem  for  the  correction  to  the  leading-order  Fourier  coefficients  may  be  written 
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down;  for  the  next-order  coefficients,  let 

u2  ~ u ^ (A9) 

and  similarly  for  the  other  variables.  Substitution  into  Eq.  (A4)  results  in 

SiwvJ1*  = ( D 2 - a2)ui1}  - \{DVou^  + DV0*u[0)),  (A10) 

2 iu(D2  - a2)41}  = ( D 2 - a2)2uP  - a2T  ^Vov^  + V*vi0)),  (All) 

2 + u{2]  = (D2  - a2)T2(0),  (A12) 

on  0 < x < 1 and  on  x0  < x < 1 we  have 

2 iufg  - L(Z)2  - a2)Tg2  = o.  (A13) 

At  the  rigid  boundaries  all  of  the  variables  are  zero,  and  at  the  interface  we  must  have 

41}  = Du{21]  = v{3l)  + DV0h{2]  + DV'h^  = 0,  (A14) 

f2(1)  + h,P  = fjg  + gshP  = (A15) 

2 iu£h{2l)  + OT2(1)  - qDT$  = 0.  (A16) 

The  term  h ^ is  identically  zero;  this  is  shown  by  a similar  procedure  to  that  for  h^2\ 
It  appears  from  this  problem  that  the  interfacial  temperature  coefficient  T2^  will  be 
nonzero  at  the  interface  and  numerical  results  show  that  the  h2  = 0(Pr-2)  and  thus  that 
h 2^  ^ 0.  These  two  observations  suggest  that  there  is  a boundary  layer  of  width  Pr1,/2  in 
the  coefficient  for  the  solid  temperature  coefficient  T^.  If  we  rescale  that  equation  with 
x = Pr1,/2x  the  equation  for  the  inner  part  of  the  solid  temperature  coefficient  is 

2iut®  - D2T = 0.  (A1Y) 

which  has  a solution  of  the  form 

T$  = Ci  exP[(!  ~ i)Vux/PT1/2}.  (A18) 

There  is  no  boundary  layer  at  the  rigid  wall,  and  this  inner  solution  matches  automat- 
ically with  the  trivial  interior  solution.  The  rest  of  the  solution  at  this  order  can  be 
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found  though  the  solution  may  no  longer  be  found  sequentially.  This  boundary  layer  is 
compatible  with  ^ 0 and  h ^ ^ 0. 

We  believe  that  the  small  or  non-existent  numerical  amplitudes  of  the  higher  harmonic 
modes  contributes  to  the  success  of  the  asymptotic  approximation. 
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Figure  Captions 

Figure  1 Schematic  diagram  of  the  crystalline  inner  annulus  (labeled  “S”)  surrounded 
by  the  liquid  phase  (labeled  “L”).  In  the  base  state  the  unperturbed  crystal-melt 
interface  is  cylindrical,  with  Ri(z',t')  = R^. 

Figure  2 Multiple-exposure  profile  of  the  instability  of  a temporally-modulated  crystal- 
melt  interface,  with  a heated  outer  cylinder  and  a cooled  inner  cylinder.  The 
material  is  succinonitrile  with  Prandtl  number  Pr  = 22.8.  An  elongated,  flexible 
particle  circulates  in  the  fluid  and  delineates  the  structure  of  the  flow  field  in  half  of 
a Taylor  vortex  cell.  The  half-cell  containing  the  particle  is  about  14  mm  from  crest 
to  trough  in  the  axial  direction,  corresponding  to  a wavenumber  of  about  a = 1.1  . 

Figure  3 Neutral  curves  in  the  (Ta,  a)-plane  for  several  values  of  the  Prandtl  number 
are  displayed.  The  other  parameter  values  are  u = 28.9,  C — 409,  q — 1.009,  Pr  = 
Ps  = 22.8,  7 = 6.2  x 10~6,  rj  = 0.69,  rjs  — 0.286.  The  uppermost  curve  ( Pr  = 0.1) 
approaches  the  rigid- walled  neutral  curve  displayed  in  [25].  The  base  state  described 
in  Section  3.1.1  is  stable  below  the  curves  and  unstable  above  them. 

Figure  4 Neutral  curves  in  the  (T,a)-plane,  where  T = PrTa.  The  parameter  values 
are  the  same  as  for  Figure  3.  The  case  Pr  = oo  is  the  result  of  the  asymptotic 
analysis  described  in  Section  6.1.1.  The  base  state  described  in  Section  3.1.1  is 
stable  below  the  curves  and  unstable  above  them. 

Figure  5 Neutral  curves  in  the  (T,a)-plane,  where  T = PrTa.  The  parameter  values 
are  the  same  as  for  Figure  3,  except  that  rj  is  no  longer  a parameter  in  the  problem. 
The  case  Pr  = oo  is  the  result  of  the  asymptotic  analysis  described  in  Section  6.1. 
The  base  state  described  in  Section  3.2.1  is  stable  below  the  curves  and  unstable 
above  them. 

Figure  6 Neutral  curves  in  the  (T,a)-plane,  where  T = PrTa,  for  several  gap  widths  77. 
All  other  parameters  are  the  same  as  Figure  4,  and  Pr  = Ps  = 22.8.  The  lowermost 
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curve  is  the  narrow  gap  limit;  the  other  curves  are  wide  gap  results. 

Figure  7 Neutral  curves  in  the  (T,u;)-plane,  where  T = PrTa,  displaying  the  wide  gap 
results  in  Table  3.  The  upper  curve  is  the  result  of  Floquet  theory  data  while  the 
lower  curve  is  the  result  of  the  asymptotic  theory.  The  base  state  described  in 
Section  3.1.1  is  unstable  above  the  curves  and  stable  below  them. 
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